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We conduct numerical experiments to investigate the spatial clustering of particles and 
bubbles in simulations of homogeneous and isotropic turbulence. Varying the Stokes pa- 
rameter and the densities, striking differences in the clustering of the particles can be 
observed. To quantify these visual findings we use the Kaplan- Yorke dimension. This local 
scaling analysis shows a dimension of approximately 1.4 for the light bubble distribution, 
whereas the distribution of very heavy particles shows a dimension of approximately 2.4. 
However, clearly separate parameter combinations yield the same dimensions. To over- 
come this degeneracy and to further develop the understanding of clustering, we perform 
a morphological (geometrical and topological) analysis of the particle distribution. For 
such an analysis, Minkowski functionals have been successfully employed in cosmology, 
in order to quantify the global geometry and topology of the large-scale distribution of 
galaxies. In the context of dispersed multiphase flow, these Minkowski functionals - be- 
ing morphological order parameters - allow us to discern the filamentary structure of the 
light particle distribution from the wall-like distribution of heavy particles around empty 
interconnected timnels. 



1. Introduction 

Even in homogeneous and isotropic turbulence particles, drops, and bubbles (all from 
now on called "particles" ) do not distribute homogeneously, but cluster, sec Crowe et al. 
(1996) for a classical review article. The clustering has strong bearing on so diverse issues 
such as aerosols and cloud formation (Falkovich et al. (2002); Celani et al. (2005); Vail- 
lancourt et al. (2002)), plankton distribution in the deep ocean (Malkiel et al. (2006)), 
sedimentation and CO2 deposition in water (Upstill-Goddard (2006)). Considerable ad- 
vances in particle tracking velocimetry (Porta et al. (2001); Hoyer et al. (2005); Bourgoin 
et al. (2006); Bewley et al. (2006); Ayyalasomayajula et al. (2006)) and in numerics (El- 
ghobashi & Truesdell (1992, 1993); Wang & Maxey (1993); Boivin et al. (1998); Druzhinin 
& Elghobashi (2001); Marchioh & Soldati (2002); Mazzitelh et al. (2003a,&); Biferale et al. 
(2005); Chun et al. (2005); van den Berg et al. (2006); Bee et al. (20066, o)) now allow 
for the acquisition of huge data sets of particle positions and velocities in turbulence. 




Figure 1. (color) Snapshots of the particle distribution in the turbulent flow fleld for St = 0.6 
for (a) /3 = 3 (bubbles), (b) (3 = 1 (tracers), and (c) (3 = (heavy particles), all for Re\ = 75. 
Corresponding videos are shown in the accompanying material (or can be downloaded from 
http: / / cfd.cineca.it /cfd /gallery / movies / particles- and-bubbles) . 



In our numerical experiments the fluid flow is simulated by the incompressible Navier- 
Stokes equation on a 128'^, a 512'^, and a 2048'^ grid with periodic boundary conditions 
and a large scale forcing, achieving Taylor- Reynolds numbers of Re\ = 75, 180 and 
400, respectively. One-way coupled point-particles are included which experience inertia 
forces, added mass forces, and drag, i.e., the particle acceleration is given by (Maxey & 
Riley (1983); Gatignol (1983)) 

^ = P§-uixit),t) - ^^{v - u{xit),t)), (1.1) 

where v = dx/dt is the particle velocity and u{x{t),t) the velocity field. Equation (1.1) 
holds in the limit of small (^ 1) particle Reynolds number and for particles whose size is 
small as compared to the Kolmogorov length scale r] - for finite size particles the results 
will naturally differ. Lift, buoyancy, two-way-coupling and particle-particle interactions 
are ignored. Next to the Reynolds number the control parameters are the particle radius 
a, the density of the fluid pf and of the particle pp. The dimensionless numbers used to 
characterize the particle in the turbulent How are the density ratio P = 3pf /{pf + 2pp) 
and the Stokes number St = Tp/r^ = o? /[ZPr]^), where Tp = a? /{2>(5v) is the particle time 
scale, V the viscosity, and r,, the Kolmogorov time scale. In our numerical study we treat 
/3 and St as independent parameters. The case /3 = (and finite St) corresponds to very 
heavy particles and /3 = 3 to very light particles, e.g. (coated) bubbles (i.e., with no-slip 
boundary conditions at the interface). (3=1 means neutral tracers and if in addition 
5*^ = we have fluid particles. 

In figure 1 snapshots of the resulting particle distributions from a simulation with 
Rex = 75, St = 0.6 and the total number of particles, N = 10^ are shown. For /3 = 
(heavy particles) and /3 = 3 (light bubbles) we observe clustering - but of different type. 
No clustering is observed for the neutral tracers with (3=1. 



2. Dimensionalities of particle and bubble clusters 

How to quantify and characterize the clustering? One way borrowed from dynamical 
system theory is to calculate the Kaplan- Yorke dimension Dky in the six-dimensional 
space spanned by the particle positions and their velocities, see Bee (2003, 2005). The 
Kaplan- Yorke dimension follows from the Lyapunov exponents (Eckmann & Ruelle 
(1985)), \i [i = 1,. . . ,6), and quantifies how contracting the dynamical system is. It is 
defined as Dy^y = K -{-Y^f^i '^il\^K+i\, K being the largest integer such that J2^i ^ 0- 
Details on the method can be found in Bee (2005) , where the case (3 = was studied. 
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Figure 2. (color) The Kaplan- Yorke dimension Dry and the correlation dimension D2 (bottom) 
of the particle distribution as function of St and {3. The black line in Dky marks the value 
Dky{P, St) = 3, the red line the value 2. For the Dky plot 10^ particles were integrated along 
the tangent space and for the D2 plot 5 • 10^ particles were followed. In both cases the averaging 
time was tens of large eddy turnovers. The particles were grouped in about 500 different types 
characterized by their (/3, St)-values. Rex = 75. 

In fig. 2, upper, the full landscape of Dky as function of /3 and St is shown for 
Re\ = 75; in fig. 3 cuts through this landscape for fixed f3 = 0, 1, and 3 are shown for 
both Re\ ~ 75 and Rex ~ 180, revealing at most a minute Reynolds number dependence 
of Dky- 

We now discuss the dependence Dky{P: St): Point-particles {St = 0, /3 = 1) do not 
experience any drag {v = u) and accordingly Dky = 3. For fixed (3 the contraction 
is strongest for a Stokes parameter St in the range between 0.5 and 1.5, i.e., when the 
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Figure 3. The Kaplan- Yorke dimension Dky vs. St for three [3 values; (3 = (o), 1 (A), 3 (□). 
Results at two different Rex are reported: Rex — 75 (filled symbols), Rex = 180 (open symbols). 
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Figure 4. (color) Projection of the particles in a slice of 19?; thickness and side 640»7 x 64077 with 
Rex = 180, for (a) St = 0.1, (b) St = 0.6, and (c) St = 4. Heavy particles {(3 = 0) are shown in 
red, bubbles (/3 = 3) in blue. To illustrate the different clustering, particles and bubbles in the 
same simulation are plotted on top of each other - note that they do not mutually interact, i.e., 
collisions and feedback on the fluid are neglected. 



characteristic time scale Tp of the particle roughly agrees with the Kolmogorov time 
scale T^. For smaller St the particles follow the small scale turbulent fluctuations and 
the clustering decreases and Dky increases. For large St the particles are so big that 
they average out the small scale turbulent fluctuations (see figure 4); correspondingly, the 
clustering decreases and Dky increases. The dynamical evolution of the heavy particles 
{P = 0) shows the strongest contraction for St « 0.5 with Dky = 2.6. Also that of the 
light particles (/3 = 3) has the strongest contraction for St ?» 0.5. Now even Dky = 1-4 
is achieved. Figures 2 and in 3 reveal another feature of Dky which is worth being 
mentioned: Dky (St, (3 = 1) slightly increase with increasing St, reflecting the fact that 
large particle, even if neutrally buoyant, do not exactly follow the flow, see the discussion 
in Barbiano et al. (2000). 

An analysis of the spatial distribution performed with the correlations dimension D2 
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(Grassberger & Procaccia (1984)) yields very similar results as that for Dky, see fig. 

2, lower (see also Bee et al. (2007) for further details on this observable). The D2 is 
computed by fitting the separation probability P2{r), whieh is the probability that the 
distance between two particles is less than r. It is assumed that P2{r) ~ r^^ as r — > 0. 
As the numerical convergence is more difficult to achieve for D2, the curve D2{P,St) is 
slightly less smooth as compared to Dxy{P, St). The finiteness of the particle samples, 
here roughly 20 snapshots of 10^ particles, leads to noisy measures especially for the 
highly clustered cases (/3 ~ 3): Point-particles may reach vanishingly small distances 
which makes it difficult to perform a proper histogram binning, and thus to compute 
statistically stable P2{r)- 

Due to its local character a metric measurement with the Kaplan- Yorke or correlation 
dimension cannot supply us with global morphological information. Indeed, as seen from 
figure 2, clusters of heavy and light particles can have the very same Dky, though they 
look very different. Indeed, the striking morphological differences of the particle and 
bubble distribution within the same turbulent velocity field are illustrated in figure 4. 

3. Morphology of particle and bubble clusters 

With global geometrical and topological order parameters we are able to distinguish 
the clustering on spatially extended, eventually interconnected, sheets from clouds or 
filamentary clustering. Consider the union set A,. = UiLo'^i'(^i) '^^ balls of radius r 
around the N particles at positions Xi^ i = 1, 2, . . . iV, thereby creating connections be- 
tween neighboring balls. The global morphology of the union set of these balls changes 
with radius r, which is employed as a diagnostic parameter. It seems sensible to re- 
quest that global geometrical and topological measures of e.g. Ar are additive, invariant 
under rotations and translations, and satisfy a certain continuity requirement. With 
these prerequisites Hadwiger (1957) proved that in three dimensions the four Minkowski 
functionals Vfj,{r), fi = 0,1,2,3, give a complete morphological characterization of the 
body Ar- The Minkowski functional Vo{r) simply is the volume of A,-, Vi{r) is a sixth 
of its surface area, V2{r) is its mean curvature divided by 37r, and V3(r) is its Euler 
characteristics. Volume and surface area are well known quantities. The integral mean 
curvature and the Euler characteristic are defined as surface integrals over the mean and 
the Gaussian curvature respectively. This definition is only applicable for bodies with 
smooth boundaries. In our case we have additional contribution from the intersection 
lines and intersection points of the spheres. Mecke et al. (1994) discuss the definitions 
of the integral mean curvature and the Euler characteristic for unions of convex bodies. 
The Euler characteristic as a topological invariant allows for several other definitions, 
like X = #(isolated bodies) — ^^(tunnels) -I- ^(completely enclosed cavities). Minkowski 
functionals have been developed for the morphological characterization of the large scale 
distribution of galaxies by Mecke et al. (1994) and have successfully been used in cos- 
mology by Kerscher et al. (1997); Kerscher (2000); Kerscher et al. (2001), porous and 
disordered media by Arns et al. (2002), dewetting phenomena by Herminghaus et al. 
(1998), statistical physics in general (Mecke (2000)) and have been recently employed to 
study magnetic structure in small scale dynamos (Wilkin et al. (2007)). 

There is an interesting relation to fractal dimensions. The scaling behaviour of the 
volume density for r — > gives the Minkowski-Bouligand dimension of a Fractal which 
equals the box-counting dimension, which is again an estimate of the Hausdorff dimension 
(see e.g. Falconer (1990)). However due to discreteness effects the scaling behaviour of 
the volume only yields a noisy estimate of the dimension. The D2 and Dky, as used 
in our case, are much more reliable measures for the dimension. Since the Minkowski 
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Figure 5. (color) Volume densities of the Minkowski functionals Vfj,{r), fj, = 0,1,2,3, for 
passive tracers (/? = 1, black dotted, corresponding to those of a Poisson distribution Mecke & 
Wagner (1991)), heavy particles (/3 = 0, red solid), and bubbles (/3 = 3, blue dashed). In all 
cases St = 0.6 and Rex = 78. Distances are in Kolmogorov units (tj). To estimated the errors 
in the Minkowski functionals we looked at several realizations of a Poisson process. The one 
sigma error bars are smaller than the linethickness of the shown curves which illustrates the 
robustness of the Minkowski functionals. - The code used for the calculations of the Minkowski 
functionals is an updated version of the code developed by Kerscher et al. (1997), based on the 
methods outlined in Mecke et al. (1994). The code is made available to the general public via 
http:/ /www. mathematik.uni-muenchen.de/~kerscher/software/ . 



functionals are dimensional quantities, one is able to define other scaling dimensions 
beyond the volume. This has been formalized and investigated for special models by 
Mecke (2000). It remains an open question how to apply this to points decorated with 
spheres. 

In figure 5 we show the volume densities of the four Minkowski functionals Wp(r) = 
V^{r)/L^, fi = 0,1,2,3, determined from the particle distributions with St = 0.6 for 
the three cases with /3 = 3 ("bubbles"), (3 = (heavy particles), and /3 = 1 (neutral 
tracers). For the neutral tracers {(3=1) the functionals coincide with the analytically 
known values (Mecke & Wagner (1991)) for randomly distributed objects. As the radius 
increases, the volume is filled until reaching complete coverage where the volume density 
Vo{r), i.e. the filling factor, reaches unity. This increase is considerably delayed for heavy 
particles and even more for bubbles, which is a characteristic feature of a clustering 
distribution produced from the empty space in between the clusters. 

The density of the surface area, measured by wi(r), increases with the radius r. As 
more and more balls overlap the growth of vi (r) slows down and the surface area reaches 
a maximum. For large radii the balls fill up the volume and no free surface area is left. 
For both the bubbles and the heavy particles the maximum of vi (r) is smaller compared 
to the Poisson case. For intermediate and large radii we observe the skewed shape of 
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Vi (r) with a significant excess of surface area on large scales compared to Poisson. The 
particles cluster on clumpy, filamentary and sheet like structures and the surface area of 
the balls is growing into the empty space in between. Especially for the bubbles (/3 — 3) 
the maximum of vi (r) is attained for considerably larger r compared to the distributions 
with /? = 0, 1, suggesting mainly separated filamentary shaped clusters. 

The density of the integral mean curvature V2 {r) allows us to differentiate convex from 
concave situations. For small radii the balls are growing outward. The main contributions 
to the integral mean curvature is positive, coming from the convex parts. Increasing the 
radius further we observe a maximum for all three cases, but as with the surface area, the 
amplitude of the maximum is reduced for heavy particles and especially for the bubbles. 
For tracer particles (/3 = 1) the empty holes start to fill up and the structure is growing 
into the cavities. Now the main contribution to the integral mean curvature is negative, 
stemming from the holes and tunnels through the structure. This concaveness is less 
pronounced for the heavy particle case {(3 = 0). Typically interconnected networks of 
tunnels show such a reduced negative contribution. For bubbles {(3 ~ 3) Ar is hardly 
concave, i.e., no holes and no tunnels develop, just as expected from isolated (filamentary) 
clusters. For large radii the balls fill up the volume, no free surface and hence no curvature 
is left. 

The topology undergoes a number of changes which we measure with the Euler char- 
acteristic. For small radii r w the balls remain separated and the volume density of 
the Euler characteristic V3{r « 0) equals the number density of the particles. As the ra- 
dius increases, balls join and the Euler characteristic decreases. Both for heavy particles 
(/3 = 0) and especially for bubbles (/? = 3) the decrease of W3(r) with increasing r is more 
dramatic, due to the clustering. When further increasing the radius r, more and more 
tunnels start to form resulting in a negative 1)3 (r). This is observed for neutral particles 
{P = 1) and less pronounced for the heavy particles (/3 = 0). No tunnels seem to form 
in the bubble distribution {/3 = 3). For neutral particles this behavior reaches a turn- 
ing point when these tunnels are blocked to form closed cavities and a second positive 
maximum of W3(r) can be seen. Neither bubbles nor heavy particles show a significant 
positive V3{r). 

To develop a physical picture we first look at the heavy particles. One expects that the 
particles are dragged out of the turbulent vortices gathering on structures around the 
vortices. A dimension Dky ~ 2.4 is indicating "fat" locally two dimensional structures. 
From the negative V2{r) for large r we conclude that these 2d-structures are mainly 
concave, i.e. growing inward. Hence these 2d-structures are not small 2d-patches, but 
enclose the vortices. The particles around the vortices form the enclosure of tunnels but 
they do not fully enclose the vortices (they do not form the "casing of a sausage"). 
Completely anclosed cavities would leed to a V3{r) larger zero for large r, which we do 
not observe. There have to be regions connecting the vortices devoid of any particles. 
In other words there exists an empty interconnected network of tunnels surrounded by 
particles. 

For the bubble distribution one expects that the bubbles are sucked into the vortices. 
The dimensional analysis with Dky ~ 1.4 is now suggesting "fat" locally one dimensional 
structures. For the bubble distribution V2{r) is almost everywhere positive, hence the 
Id-structures remain convex for all radii. Similar V2{r) is positiv on all scales. Hence, 
neither tunnels devoid of bubbles nor empty cavities enclosed by bubbles form. The Id- 
structures built from the bubbles are separated filamentary clusters. This also shows that 
the bubbles do not form a percolating structure throughout the simulation. 

Minkowski functionals calculated from points decorated with spheres do depend on the 
number density of the point distribution. One can derive explicit expression in terms of 



8 



Calzavarini, Kerscher, Lohse, and Toschi 




r [rj] r [77] 

Figure 6. (color) Volume densities of the Minkowski functionals v^{r), = 0, 1, 2, 3 for the data 
set with /3 — and St = 0.6. We compare the Minkowski functionals for samples with a different 
number of points, obtained by randomly thinning the original data set with 100k points. 100k 
points (black solid), 90k points (red solid), 80k points (blue solid), 50k points (black dotted), 
10k points (red dotted), Ik points (blue dotted). 

high order correlation functions quantifying the non-trivial dependence on the number 
density (sec e.g. Mecke (2000)). In our analysis we always compare data sets with the same 
number of points. However one may ask the question how well our Minkowski functionals 
are converged. In Fig. 6 the Minkowski functionals for a randomly subsampled data set 
are shown. Considering only 80% of the points we obtain nearly identical result as for the 
full sample. This can be understood as follows: We decorate the particles with spheres. 
In our case the particles follow well defined structures. For a radius equal to zero, V3 
equals the number density. That is why the U3 curves fan out for small radii. Already for 
r » 0.577 the curves for 100%, 90% and 80% overlap again. In this sense our results are 
well converged. 

Obviously, the convergence must get lost if we further and further reduce the number of 
points. Randomly subsampling (thinning) will inevitably lead to a convergence towards 
Poisson (see the blue dotted line in fig. 6) - the discreteness effects become more and 
more important. 

Up to now we investigated three extreme cases with St = 0.6 and (3 — 0, 1,3. These 
three particle distributions also show different Dky , hence their morphological differences 
do not come as a surprise. Now we investigate the morphology of the particle distributions 
with similar Dxy but quite different physical parameters /3 and St. In figure 7 we compare 
heavy particles (/3 = 0, St = 0.6) with Dky = 2.59 with Hght bubbles (/3 = 3, St = 0.103) 
with Dky = 2.55. Clearly, their morphological properties differ. Although both have 
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Figure 7. (color) Volume densities of Minkowski functionals ^^(r), ^ = 0, 1, 2, 3 for heavy par- 
ticles with /? = and St = 0.6 (red solid) and for light particles with /3 = 3 and St = 0.103 (red 
short dashed): Whereas for these two particle types Dky basically is the same, the Minkowski 
functionals clearly differ. As a second pair with nearly identical Dky but different Minkowski 
functionals we compare two light particle distribution with differing Stokes parameter: /3 = 3, 
5** = 0.4 (green long dashed) against {3 = ?>, St = 1.75 (green dashed-dotted) . Poisson distribu- 
tion behavior is reported for reference (black dotted). 

similar Dky, the distribution of the heavy particle shows a stronger deviation from the 
Poisson distributed points than the light bubbles. As another example, in figure 7 we 
also compare the distribution of two sets of light bubbles, both with /3 = 3 and identical 
Dky — 1.65, but with different Stokes number. The data set with St = 1.75 shows a 
distribution dominated by isolated clusters whereas the data set with St = 0.4 allows for 
some interconnected empty tunnels. 

We now address the Reynolds number dependence of the results from the morpholog- 
ical analysis. We compare the results for particles in turbulent flow with Re\ = 75 with 
simulations at higher Reynolds numbers Re\ = 180 and Re\ = 400, keeping the param- 
eters (3 = and St = 0.6 fixed. Only a weak dependence of the Minkowski functionals on 
the turbulence intensity can be seen in Fig. 8. A higher Reynolds number consistently 
corresponds to a slightly more pronounced deviation from the Poisson behavior. We ob- 
serve similar results for the bubble distribution with (3 = 3, St = 0.6. This result agrees 
with the very weak Reynolds dependence of Dky already revealed in Fig. 3 and in Bee 
et al. (2006a). 

4. Conclusions and outlook 



Often the small-scale behavior in a turbulent dispersed multiphase flow is of interest. 
Then scaling or contraction indices like the Kaplan- Yorke Dimension Dky are the main 
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Figure 8. (color) Volume densities of the Minkowski functionals v^{r), = 0, 1, 2, 3 for particles 
with /3 = and St = 0.6 in simulations with different Reynold numbers: Re\ = 75 (red solid), 
Rex = 180 (blue short dashed) , Rex — 400 (blue long dashed) . Poisson distribution behavior is 
reported for reference (black dotted). All point sets have been subsampled to the same number 
density. 



tools. With Dky as a function of (3 and St a simple picture emerges: For heavy particles 
(/? < 1) we observe in the extreme cases a scaling reminiscent of local planar structures, 
for hght bubbles (/3 > 1) we find indications for local linear structures. The strength of 
the deviation from Dky = 3 is further modulated by the Stokes parameter St. If one is 
interested in the connectivity and other (global) morphological features of the particle 
distribution, the picture becomes significantly more complex. The distribution of heavy 
particles seems to allow interconnected empty tunnels, whereas the distribution of bub- 
bles typically shows isolated filamentary structures. But, depending on both f3 and St, 
also interconnected structures appear in the bubble distribution. Using Minkowski func- 
tionals as morphological order parameters we are able to quantify these geometrical and 
topological features in a unique way. Especially, we are able to overcome the degeneracy 
seen in Dky as a function of /3 and 5*^. 

The next step is to use these tools to compare data obtained by parti- 
cle tracking velocimetry (PTV) with those obtained from numerical simulations, 
in order to test whether effective force models lead to the same particle dis- 
tribution as measured in experiment. To facilitate such an analysis, we make 
the code to calculate Minkowski functionals available to the general public, see 
http://www.mathematik.uni-muenchen.de/~kerscher/software/. Raw data can be down- 
loaded from the International CFD database, iCFDdatabase (http://cfd.cineca.it). 
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